////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                //
//  Ebrahim Foulaadvand 24 August 2011                             //
//                                                                                                //
//  This code gives the few first eigen values/eigenvectors for the time-independent Schroedinger //
//  equation for an infinite square well by the shooting method.                                  //
//                                                                                                //
//                                                                                                //
////////////////////////////////////////////////////////////////////////////////////////////////////


#pragma hdrstop
#include <iostream>
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <values.h>
#include <vector.h>
#include <list.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

using namespace std;

main(){

int i,N=100,exit;

double m=1.0,hbar=1.0,a=1.,E=4.9,delx,cutoff,delE=0.001,precision=0.01,terminatedphi,phiRboundary;

ofstream file1 ("V.plt");
ofstream file2 ("PsiR E=40.93.plt");
ofstream file3 ("PsiL E=40.93.plt");
ofstream file4 ("PsiR E=6.55.plt");
ofstream file5 ("PsiL E=6.55.plt");
ofstream file6 ("PsiR E=6.65.plt");
ofstream file7 ("PsiL E=6.65.plt");
ofstream file8 ("Ground stateL.plt");
ofstream file9 ("Ground stateR.plt");


delx=a/double(N); cout<<"delx=  "<<delx<<endl;

vector<double> V(N+1,0),phi(N+1,0);


//phi[0]=1; // for even eigenstates (zero derivative at x=0) therefore phi[-1]=phi[0];
phi[0]=0; // We consider odd eigenstates

phi[1]=1*delx; // for odd states


//cutoff=5*phi[0];  //even states
cutoff=5*phi[1];  //odd states


for (int i=0;i<=N;i++) { V[i]=0.; }


phiRboundary=2*precision;

while (fabs(phiRboundary)>precision){

terminatedphi=phi[0]; exit=0; i=0;

//phi[1]=2*phi[0]- phi[0] + 2*m*(delx/hbar)*phi[0]*(delx/hbar)*(V[0]-E); // for even states

phi[1]=1*delx; // for odd states

while (i<N-1){ i++; //cout<<"i= "<<i<<endl;

phi[i+1]=2*phi[i]- phi[i-1] + 2*m*(delx/hbar)*phi[i]*(delx/hbar)*(V[i]-E); //cout<<"phi["<<i+1<<"]= "<<phi[i+1]<<endl;

terminatedphi=phi[i+1];

if (fabs(terminatedphi)>cutoff) { i=2*N; exit=1; }


} // end of second while


if (exit==1) { //cout<< "terminatedphi= "<<terminatedphi<<endl;

if (terminatedphi>0)  E+=delE;  if (terminatedphi<0)  E-=delE; }

if (exit==0) { phiRboundary=phi[N]; //cout<<"phRboundary= "<<phiRboundary<<endl; getchar();

if (phiRboundary>0)  E+=delE; if (phiRboundary<0)  E-=delE;

}

//cout<<"E= "<<E<<endl; getchar();


} // end of first while


 for (int i=0;i<=N;i++) { file8<<-(i*delx)<<"  "<<phi[i]<<endl; file9<<+(i*delx)<<"  "<<phi[i]<<endl; }

 cout<<"Ground state energy= "<<E<<endl; getchar();


/*
E=9.10; cout<<"E= "<<E<<endl; getchar();

phi[1]=2*phi[0]- phi[0] + 2*m*(delx/hbar)*phi[0]*(delx/hbar)*(V[0]-E);  //cout<<"phi[1]= "<<phi[1]<<endl;

file2<<(0*delx)<<"  "<<phi[0]<<endl;
file2<<+(1*delx)<<"  "<<phi[1]<<endl;
file3<<-(1*delx)<<"  "<<phi[1]<<endl;


for (int i=1;i<N;i++) {

phi[i+1]=2*phi[i]- phi[i-1] + 2*m*(delx/hbar)*phi[i]*(delx/hbar)*(V[i]-E); //cout<<"phi["<<i+1<<"]= "<<phi[i+1]<<endl; getchar();

file2<<+(i+1)*delx<<"  "<<phi[i+1]<<endl;
file3<<-(i+1)*delx<<"  "<<phi[i+1]<<endl;

}
*/





}// main